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Abstract 

We assess the validity of a single step Godunov scheme for the solution of the magneto-hydrodynamics 
equations in more than one dimension. The scheme is second-order accurate and the temporal discretization 
is based on the dimensionally unsplit Corner Transport Upwind (CTU) method of Colella. The proposed 
scheme employs a cell-centered representation of the primary fluid variables (including magnetic field) and 
conserves mass, momentum, magnetic induction and energy. A variant of the scheme, which breaks momen- 
tum and energy conservation, is also considered. Divergence errors are transported out of the domain and 
damped using the mixed hyperbolic/parabolic divergence cleaning technique by Dedner et al. (J. Comput. 
Phys., 175, 2002). The strength and accuracy of the scheme are verified by a direct comparison with the eight- 
wave formulation (also employing a cell-centered representation) and with the popular constrained transport 
method, where magnetic field components retain a staggered collocation inside the computational cell. Re- 
sults obtained from two- and three-dimensional test problems indicate that the newly proposed scheme is 
robust, accurate and competitive with recent implementations of the constrained transport method while 
being considerably easier to implement in existing hydro codes. 

Key words: Magnetohydrodynamics, Compressible Flow, Unsplit scheme, High-order Godunov method, 
Cell-centered method 



1. Introduction 

A primary aspect in building stable and robust Godunov type schemes for the numerical solution of the 
compressible magnetohydrodynamics (MHD) equations relies on an accurate way to control the solenoidal 
property of the magnetic field while preserving the conservation properties of the underlying physical laws. 
Failure to fulfill cither requisite has been reported as a potential hassle leading to unphysical effects such 
as plasma acceleration in the direction of the field, incorrect jump conditions, wrong propagation speed of 
discontinuities and odd-even decoupling, see 26|, ij]. A comprehensive body of literature has been dedicated 



to this subject and several str ateg ies to enforce the V ■ B = condition in Godunov-type codes have been 
proposed, see for example [28, 2J, 25|, 0, [2(| and, more recently, 0,0, 15, 23, Tj^. The robustness of one 
method over another can be established on a practical base by extensive numerical testing, see 0,0). 

In a first class of schemes, the magnetic field is discretized as a cell-centered quantity and the usual for- 
malism already developed for the Euler equation can be extended in a natural way. Cell-centered methods 
are appealing since the extensions to adaptive and/or unstructured grids are of straightforward implementa- 
tion. Moreover, the same interpolation scheme and stencil used for the other hydrodynamic variables can be 
easily adapted since all quantities are discretized at the same spatial location, thus facilitating the extension 
to schemes possessing higher than second order accuracy. Unfortunately, numerical methods based on a 
cell-centered discretization do not naturally preserve Gauss's law of electromagnctism, even if V ■ B = 
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initially. In the approach suggested by Powell 2l|, |22j, Gauss's law for magnetism is discarded in the deriva- 
tion of the MHD equations and the resulting system of hyperbolic laws is no longer conservative by the 
appearance of a source term proportional to V • B. Although the source term should be physically zero at 
the continuous level, Powell showed that its inclusion changes the character of the equations by introducing 
an additional eighth wave corresponding to the propagation of jumps in the component of magnetic field 
normal to a given interface. A different approach is followed in the projection scheme (5, l28l. I24I [3| . where a 
Helmholtz-Hodge decomposition is applied to resolve B as the sum of an irrotational and a solenoidal part, 
associated with scalar and vector potentials. A cleaning step allows to recover the divergence- free magnetic 
field by subtracting the unphysical contribution coming from the irrotational component at the extra cost of 
solving a Poisson equation. In the approach of Dedner et al. [IH, the divergence- free constraint is enforced 
by solving a modified system of conservation laws where the induction equation is coupled to a generalized 
Lagrange multiplier. Dedner et al. showed that the choice of mixed hyperbolic/parabolic correction offers 
both propagation and dissipation of divergence errors with the maximal admissible characteristic speed, 
independently of the fluid velocity. This approach preserves the full conservation form of the original MHD 
system at the minimal cost of introducing one additional variable in the system and will be our scheme of 
choice. Finally, Torrilhon 27 1 (see also [H) showed a general procedure to modify the inter-cell fluxes in the 
framework of a flux distribution scheme that preserves the value of a certain discrete divergence operator in 
each control volume. 

A different strategy is followed in the constrained transport (CT) methods, originally devised by [l3| 
and later built into the framework of shock-capturing Godunov methods by a number of investigators, e.g., 
@i B EH, 15, EH- In this class of schemes, the magnetic field has a staggered representation whereby the 
different components live on the face they are normal to. Hydrodynamic variables (density, velocity and 
pressure) retains their usual collocation at the cell center. CT schemes preserve the divergence-free condition 
to machine accuracy in an integral sense since the magnetic field is treated as a surface averaged quantity and 
thus more naturally updated using Stokes' theorem. This evolutionary step involves the construction of a 
line-averaged electric field along the face edges, thereby requiring some sort of reconstruction or averaging of 
the electromotive force from the face center (where different components are usually available as face centered 
upwind Godunov fluxes) to the edges. A variety of different strategies have been suggested, including simple 
arithmetic averaging |25| , solution of 2-D Riemann problems [l8 , 14 1 or other somewhat more empirical 
approaches 0, Ell Tl7j ]. The staggered collocation of magnetic and electric field variables in CT schemes 
makes their extension to adaptive grids rather arduous and costly. Besides, significant efforts have to be spent 
in order to develop schemes with spatial accuracy of order higher than second. An alternative constrained 
transport method, based on the direct solution of the magnetic potential equation (thus avoiding staggered 
grids), has been presented by [23| . 

In the present work we propose a new fully unsplit Godunov scheme for multidimensional MHD, based 
on a combination of the Corner Transport Upwind of [8j and the mixed hyperbolic/parabolic divergence 
cleaning technique of [12] (CTU-GLM). The proposed scheme has second order accuracy in both space 
and time and adopts a cell-centered spatial collocation (no staggered mesh) of all flow variables, including 
the magnetic field. The scheme is fully conservative in mass, momentum, magnetic induction and energy 
and the divergence- free constraint is enforced via a mixed hyperbolic/parabolic correction which avoids the 
computational cost associated with an elliptic cleaning deriving from a Hodge projection. A variant of the 
scheme, which introduces divergence source terms breaking the conservative properties of some equations, is 
also presented. We assess the accuracy and robustness of the scheme by a direct quantitative comparison with 
the 8- wave formulation of [22| and the recently developed constrained transport method of lfj lH • Other 



similar implementations may be found in 14j, |17| . The comparison is conveniently handled using the PLUTO 



code for computational astrophysics [201 ] where both cell-centered and staggered-mesh implementations are 
available. 

Our motivating efforts are driven by issues of simplicity, efficiency and flexibility. In this sense, the 
benefits offered by a method where all of the primary flow variables are discretized at the same spatial 
location considerably ease the extension to adaptive grids, to more complex physics and to schemes with 
higher than second order accuracy. The latter possibility will be explored in a companion paper. 
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2. The Constrained GLM-MHD Equations 



In the approach of [12|, the divergence constraint of the magnetic field (Gauss's law) is coupled to 
Faraday's equation by introducing a new scalar field function or generalized Lagrangian multiplier -0. The 
second and third Maxwell's equations are thus replaced by 
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where I? is a linear differential operator. Dedner et al. built this approach into the MHD equations and 
showed that a satisfactory explicit approximation may be obtained by choosing a mixed hyperbolic/parabolic 
correction, according to which 2?(V>) = c^ 2 dt^ + c~ 2 ip where Ch and c p are constants. Direct manipulation 
of the modified Maxwell's equations |T]) leads to the telegraph equation, 
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which implies that divergence errors are propagated to the domain boundaries at finite speed Ch and decay 
with time and distance. The constant ratio c^/c 2 , which has the dimension of inverse time, sets the damping 
rate. In the limiting case of c p — > oo, one retrieves the simple hyperbolic correction and Eq. © reduces to 
an ordinary wave equation. 

The GLM-Maxwell's equations fl} can be coupled to the equations of magnetohydrodynamics writ- 
ten in their conservative form. The resulting system is called the generalized Lagrange multiplier (GLM) 
formulation of the MHD equations (GLM-MHD) and is comprised of the following nine evolution equations: 
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where p, v, p and B are the mass density, velocity, gas pressure and magnetic field, respectively. Total 
energy E and gas pressure are related by the ideal gas law, E = p/(T — 1) + pv 2 /2 + B 2 /2, where T is the 
specific heat ratio. Notice that we have conveniently switched, using vector identities, to the divergence form 
of the induction equation, more appropriate for the cell-centered finite volume formalism. The constrained 
GLM-MHD equations © are hyperbolic and fully conservative in all flow variables with the exception of 
the unphysical scalar field tp which satisfies a non-homogeneous equation with a source term. Divergence 
errors propagate with speed Ch independently of the flow velocity, thus avoiding accumulation in presence 
of stagnation points. The presence of the source term is responsible for damping divergence errors as they 
propagate. 

Dedner et al. also considered a slightly different constrained formulation, in which the Lorentz force 
term in the MHD equations is directly derived from the GLM-Maxwell equations. In this case, the system 
([3]) is extended by an additional source term on the right hand side, namely 



Seglm = [0, -(V • B)B, 0, -B • W, 0] 1 , 
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(4) 



where the non-zero entries correspond to the momentum and energy equations. Dedncr called the system 
^ augmented with the source term (0| on its right hand side the extended GLM (EGLM) formulation 
of the MHD equations. Although the system breaks conservation of energy and momentum, it still holds 
some attractive features and we found it, in our experience, a more robust scheme in presence of strong 
discontinuity propagating through highly magnetized environments. 



3. The CTU-GLM scheme 



We now illustrate the detailed steps of our new cell-centered numerical scheme. The derivation is shown 
for the conservative GLM scheme, whereas modifications relevant to the EGLM formulation are described 



We adopt a Cartesian system of coordinates and re-write the system of equations in ([3]) as 
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where d, I = x,y,z label the different component and flux contributions in the three directions while ddi 
is the delta Kronecker symbol. The system of equations given in ((SJ) is advanced in time by solving the 
homogeneous part separately from the source term contribution, in an operator-split fashion: 



where A and S are the advection and source step operators separately described in A3. H and 



(6) 

respectively. 



3.1. Advection Step 

During the homogeneous step, we adopt a numerical discretization of §5§ based on the corner trans- 
port upwind (CTU) method of [8|. For simplicity, we will assume hereafter an equally-spaced grid with 
computational cells centered in (xi,yj,Zk) having size Ax x Ay x Az. For the sake of exposition, we omit 
the subscript (i,j,k) from cell centered quantities while keeping the half increment index notation when 
referring to the interfaces, e.g., Pj+± = Pi j+i fe- An explicit second order accurate discretization of Eqns. 
(O, based on a time-centered flux evaluation, reads 
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where U = (ft, pv, B, E, ip) is the state vector of conservative variables. The expression in square brackets 
provides a conservative discretization of the divergence operator appearing in the original conservation laws 
with F, G and H being suitable numerical approximations to the flux contributions in §5§ coming from 
the I = x 7 y,z directions, respectively In the CTU approach, numerical fluxes are computed by solving a 
Riemann problem between suitable time-centered left and right states, i.e., 
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where V 



{ft, v, B,p, ip) T is the state vector of primitive variables and 7Z(-,-) denotes the flux obtained 

by means of a Riemann solver, see W6.2\ The corner-coupled states, V"^ 2 and V™^ 1 2 _, are computed via 
a Taylor expansion consisting of an evolutionary step in the direction normal to a given interface f £!3.1.ip 
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followed by a correction step involving transverse flux gradients f i)3.1.2p . The algorithm requires a total of 
6 solution to the Ricmann problem per zone per step. 

The time increment At" is computed via the Courant-Friedrichs-Levy (CFL) condition: 



At" = C a - 
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where the maximum and minimum arc taken over all zones and Cf tX , Cf,y, c f,z are the fast magneto-sonic 
speeds in the three directions, see A3. 1.11 C a is the Courant number and, for the 6-solve CTU presented 
here, is restricted to C a < 1 in two dimensions and C a < 1/2 in three dimensions. 

3.1.1. Normal Predictors 

During the computation of the normal predictors, we take advantage of the primitive (or quasi-linear) 
form of the equations. By discarding contributions from y and z and considering the reconstruction process 
in the x direction only, one has 
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is the usual matrix of the MHD equations in primitive form plus the addition of a fifth row and a ninth 
column. The source terms Sb x and S^, are of crucial importance for the accuracy of the scheme in multi- 
dimensions @, [IB, 17 1 and take the form 
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The matrix A x of the quasi-linear form is diagonalizable with the same eigenvalues as the ordinary MHD 
equations plus two new additional entries Ch and — Ch, for a total of 9 characteristic waves: 
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are the fast magneto-sonic (c/ with the + sign), slow magneto-sonic (c s with the — sign) and Alfvcn 
velocities. The two additional modes ±c n are decoupled from the remaining ones and corresponds to waves 
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carrying jumps in B x and ip. The constant Ch gives the speed of propagation of local divergence errors and 
is chosen to be the maximum speed compatible with the time step restriction, in other words 



Ch=max.(\v x \+Cf >x ,\v y \+Cfy,\v x \+cj z ) , (15) 

Finally, the corresponding left (l fc ) and right (r k ) eigenvectors are given in Appendix [A"l 

Using the characteristic decomposition of the quasi-linear form (fTH)) . we extrapolate V(ccj,£™) from the 
cell center to the edges x i± 1 for a time increment At™/2. During this step we only consider the contribution 
of those waves traveling from the center to the given interface and discard any interaction between neighbor 
cells. The resulting construction yields the normal predictors 

1 / \ k At n \ At™ 

where only positive waves (A fe > 0, k = 1, ...,9) contribute to the left of the i + ^ interface (i, +) while only 
negative waves (A fc < 0) are considered to the right of the i — \ interface («,—)• The undivided differences 
AB X and Aip may be computed using a standard centered finite difference approximation. The jump 
contribution from the k— th characteristic field is denoted with AVf = Awfrf where rf is the corresponding 
right eigenvector and Aw' is a limited slope in the k— th characteristic variable, 

Aw* = Lim fl* ■ AV" ■ AVh) , (17) 



*+_ 

where AV™ ±i = ± (V™ ±1 — Vf), l fe is the k— th primitive left eigenvector and Lim(-, •) is a limiter function, 
e.g. 

t- tx x \ signQS-) + sign(<$ + ) . / . + \ 

Lim(6_,<5+) = mm I /3|d_|,/3|d + |, j . (18) 

Usually taking (3 = 2 gives the largest compression. However, for problems involving strong shocks, we 
found setting (3 = 1 for nonlinear fields (fast and slow shocks) and (3 = 2 for the linear fields to give a more 
robust recipe. 

3.1.2. Transverse Predictors 

Once the normal predictor states have been computed, we solve a Riemann problem at constant y— and 
z— faces to obtain the transverse fluxes, e.g., 

G* + i =ft(V* + ,V* +1) _) , H * + i = ^(V*, + ,V* +1 ._) , (19) 

where left and right states have been computed during the normal predictor stages in the y and z direction. 
The solution of the Riemann problem follows the guidelines illustrated in £13.21 where the linear sub-system 
formed by the longitudinal magnetic field component and the Lagrange multiplier is preliminary solved 
before a standard 7— wave Riemann solver is applied. Transverse flux gradients are then added to the 
normal predictors (HHJ) once they are transformed back to conservative variables. This yields the corner 
coupled states: 

„+! At ( G* +i ~G*_ l2 H* + , -H* A 

U. ± i -U i±h 2 y Ay + Az j , (20) 

where U* is obtained by converting V* to conservative variables. 

We recall that the starting point in the derivation of Eq. (|20[) may be viewed, in its simplest form, as a 
first order Taylor expansion around the cell center (xi,t n ), 

^+h^r,n,d\]-Ax dV-At^( dV-Ax AtdF?\ At (dG* 5H|\ 
U *±i ~ U i ± ~dx~~ ~df^2 ~ y i ~dx~^2 2~dx )~~2 \~dy~ ~dz~ J ' (21) 

6 



where the temporal derivative dXJ/dt has been replaced, in the second expression, by taking advantage of 
the original conservation law and the different terms have been grouped according to the step in which they 
are computed (i.e., Eg 1161 and Eg I20p . In this perspective, the input states entering in the computation 
of the transverse fluxes (fl"9|) may be slightly modified by 0(At 2 ) in order to more accurately represent 
the V • B term in the construction of the scalar multiplier ip. To better understand this minor correction, 
we rewrite the ifi component of the interface states (f2"U)) in 2D using, for the sake of simplicity, a simple 
MUSCL-Hancock step during the normal predictor: 
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Clearly, the multidimensional terms approximating V ■ B in the square bracket of Eq. 
normal (AS™) and a transverse (B* x — B* 



split into a 

directional contribution. Since the first one is taken at 



time level n while the second term comes from solving a Riemann problem between normal predictors in 
the y direction (extrapolated a t n + At n /2), these contributions are not taken at the same time level but 
are spaced by At n /2. In practice, from the tests included here and several others we found evidence that a 
better balance is achieved if one replaces, in the input states of (Q33), the longitudinal field component with 
its interpolated value at time level n, i.e., B* J : ± — > _B™ ± AB™/2 (or, cquivalently with the value obtained 
by setting At = in Eq. ITH|) . Note that this is a second-order correction that does not alter the accuracy 
of the scheme and only affects the solution of the Riemann problem in computing the transverse fluxes p9[) 
but does not concern the definitions of the normal predictors. Although this is not an essential step, it was 
found to improve the accuracy in the numerical tests presented in S]4] 



3.2. Solving the Riemann Problem 

In the case of the GLM-MHD equations, left and right input states to the Riemann solver bring a 

set of 9 jumps propagating along the 7 standard characteristic MHD waves (i.e. fast, slow, rotational pairs 
and one entropy modes) as well as 2 additional modes carrying jumps only in the normal (longitudinal) 
component of B and ip. Nonetheless, when solving a one-dimensional Riemann problem at a zone interface 
(say the x direction), these additional waves are decoupled from the remaining ones and are described by 
the 2x2 linear hyperbolic system 

dB x _d^_ 
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For a generic pair of left and right input states (-B^LjV'l) and (B Xi r,iPr), the Godunov flux of the system 
(|2"3")l can be computed exactly as 
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This allows to carry out the solution of the 2x2 linear Riemann problem separately before using any 
standard 7- wave Riemann solver for the one-dimensional MHD equations. The longitudinal component of 
the magnetic field £?*, preliminary computed with (|24j) . enters hence the ordinary Riemann flux computation 
as a constant parameter. 

In other words, given the arbitrary left and right states "V l and Vr, input to the Riemann problem, we 
compute 

K(V L ,V R )=K 7 (VhV* R ) (25) 

where VJ (S = L, R) is the same as Vs with (B x ,s, ips) replaced by (B*,ip*) and H7 is a standard 7— wave 
Riemann solver. In this work, we will employ the linearized Riemann solver of Roe, in the version of Q. 
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3.3. Source Step 

During the source step we solve the initial value problem given by the last of equations © without the 
V • B term, that is, 



supplemented with the initial condition ip(°> given by the output of the most recent step. The constant c 2 
has the dimension of length squared over time and thus can be regarded as a diffusion coefficient. Dedner 
et al. prescribe an optimal value c^/ch = 0.18 independently of the mesh spacing; however, we suspect 
this definition to be incomplete, since c^/ch has the dimension of length and thus it is not a dimensionless 
quantity. Our numerical experiments indicate that divergence errors are minimized when the parameter 
a = Ahch/Cp (where Ah = min(Ax, Ay, Az)) lies in the range a £ [0,1], depending on the particular 
problem. In first approximation this value can be regarded as grid-independent although we have verified 
a weak tendency to decrease as the mesh thickens. Using the definition of a, Eq. (f2"o| can be integrated 
exactly for a time increment At 71 , yielding 

^")^(o) exp |_ a _|_j ! with a = Ah°±. (27) 

Note that, when is chosen using Eq. (|15p . the argument of the exponential becomes simply (—C a a). 
Finally, we comment out that the dimensionless a parameter can be regarded as the ratio of the diffusive 
and advective time scales, i.e., a = Atd/At a , where Atd = Ah? jc^ and At a = Ah/ch- 

3.4- Modifications for the Extended GLM (EGLM) formulation 

The extended GLM-MHD (EGLM-MHD) equations may be derived from the primitive MHD equations 
rather than the conservative ones, (T^]. In this approach, the divergence part of the Lorentz force is added 
to the momentum flux and an additional source term, given by ([2|), is introduced into the system. The 
construction of the normal predictor states carried out in ^3. 1.11 remains the same with the exception of the 
source terms (|12p which must be replaced by 

Sb. = [0,0, 0,0,0, v y ,v z ,-{T- l)v-B,0] T , = . (28) 

Since the corner coupled states in Eq. (|20[) are obtained in conservative variables, they must also be 
augmented with the source term contribution (Eq. 3]) and thus replaced by 

U i= ^"i 2 — ► + — (S EGLAIy + ^EGLM,z) ■ (29) 

Likewise, the final update Eq. (J7J) becomes 

U n+1 - U-+ 1 + At (s n E + J LM x + Sl + J LM y + S n E + J LM z ) . (30) 

In Eq. (|29p and (|30|) we have split the source term into contributions coming from the derivatives in the x, 
y and z directions. For each term we take advantage of the upwind fluxes computed in the corresponding 
direction during the Riemann solver step. For example, during the y— sweep we compute the momentum 
and energy sources in Seglm,u as 

.1>B, v ,( B l.,H- B l.i-i\ „ (*i*i-*>-! S 



where B* and ip* follows from the solution of the linear 2x2 Riemann problem The cell-centered 

magnetic field is evaluated at t n for the computation of the corner coupled states ([2"5]) and by averaging to 
cell-center the final interface values for the final update, Eq. pO]) . 



4. Numerical Tests 



We now proceed to a direct verification of the CTU-GLM and CTU-EGLM algorithms developed in the 
previous sections. A test suite of standard two- and three-dimensional MHD problems has been selected in 
order to monitor and quantify the accuracy of the proposed schemes. For the sake of comparison, we extend 
the verification process to other two well known methods, namely, Powell's eight wave formulation 
based on a cell-centered approach and the constrained transport (CT) scheme of lBL 16j| using a staggered 
formulation. The four selected algorithms, "GLM" , "EGLM" , "8W" and "CT" , have been built into the CTU 
methodology and have been implemented in the current distribution of the PLUTO code for astrophysical 
gas-dynamics (20| available at http://plutocode.to.astro.it. Adopting the same numerical framework provides 
a practical way for a convenient and extensive inter-scheme comparison. 

In the following test problems the scalar field function ip will be always initialized to zero and thus 
omitted from the definition of the initial conditions. Moreover, unless otherwise stated, the specific heat 
ratio will be set to T = 5/3 and the default Courant number is set to C a = 0.8 in two dimensions and 
C a = 0.4 in three dimensions. Errors for any flow quantity Q are computed using the L\ discrete norm 
defined by 

(Q) = tt^ftf £ I Q«,* - I ( 32 ) 



N x N y N z 



i,j,k 



where N x , N y and N z are the number of points in the three directions, Ql°j k 
summation extends to all grid zones. 



is a reference solution and the 



4-1. Propagation of Circularly polarized Alfven Waves 

Circularly polarized Alfven waves are an exact nonlinear solution of the compressible MHD equations 
thus providing an excellent code benchmark. For a planar wave propagating along the x direction with 
angular frequency u> and wave number k, the transverse components of velocity and magnetic fields trace 
circles in the yz plane and the solution can be written as 
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where <f) = kx — ut, tojk = Vq x ± c a is the corresponding phase velocity (c a = 1 is the Alfven speed) and 
A = 1/10 is the wave amplitude. The plus or minus sign corresponds to right or left propagating waves, 
respectively. The constants Vq x , Vq v , vq z give the translational velocity components in the three directions. 
Density and pressure remain constant and equal to their initial values po = 1 and po = 0.1 since torsional 
Alfven waves do not involve any compression. 

Here we consider a rotated version of the one-dimensional solution given by (|33[) and specify the orien- 
tation of the wave vector k = (k x ,k y ,k z ) in a three dimensional space x,y,z through the angles a and (3 
such that 

k k 

tan a=-r L , tan/3 = . (34) 

The full 3D solution is then recovered by rotating the original one dimensional frame by an angle 7 = 
tan _1 (cosatan/3) around the y axis and subsequently by an angle a around the z axis. The resulting 
transformation leaves scalar quantities invariant and produce vectors rotation q — > R 7Q q, where 



' cos a cos 7 — sin a — cos a sin 7 \ 

sin a cos 7 cos a — sin a sin 7 
1 sin 7 cos 7 J 



R 7 a 



/ cos a cos 7 sin a cos 7 sin 7 \ 

— sin a cos a 

V — cos a sin 7 — sin a sin 7 cos 7 J 



, (35) 
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are the rotation matrix and its inverse whereas q is a three-dimensional vector. Note that <fi is now given 
by cf> = k • x — ojt where to = |k|(t;ox i c a ) and 



|k| =fc K yl + tan 2 a + tan 2 /3 (36) 

is the wavenumber corresponding to a wavelength A = 27r/|k| and period T = 2tt/lu. 

In order to ensure correct periodicity we assume, without loss of generality, k x = 2tt and pattern the 
computational domain such that one wave period is prescribed in each grid direction, i.e., x € [0, 1], y € 
[0, 1/ tana] and z E [0, 1/ tan/3]. Also, for the tests discussed here, we consider standing waves and thus 
set vq x = vq v = vq z = 0. With these definitions the wave returns into the original position after one period 
T = X/c a with 

T = 1 . (37) 

V 1 + tan 2 a + tan 2 (3 

4-1.1. Two- Dimensional Propagation 



Table 1: Errors (in L\ norm) and orders of accuracy for the two and three-dimensional circularly polarized Alfven wave tests. 
The first and second columns refer to the numerical scheme and the number of points in the x direction. Columns 3-4 and 
5-6 show the result obtained in the 2D problem with Courant number of C a = 0.8 and C a = 0.4, respectively. The last two 
columns corresponds to the three dimensional case. 







2D, C a 


= 0.8 




2D, C a 


= 0.4 


3D, C a 


= 0.4 


Scheme 


N x 


L\ Error 


Li 


order 


L\ Error 


L\ order 


L\ Error 


L\ order 


GLM 


16 


2.46E-002 






2.60E-002 




3.19E-002 






32 


4.56E-003 




2.43 


5.17E-003 


2.33 


5.66E-003 


2.50 




64 


1.16E-003 




1.97 


1.27E-003 


2.03 


1.15E-003 


2.30 




128 


3.19E-004 




1.87 


3.02E-004 


2.07 


3.03E-004 


1.92 




256 


8.48E-005 




1.91 


7.01E-005 


2.11 


8.05E-005 


1.91 


CT 


16 


2.54E-002 






2.79E-002 




3.44E-002 






32 


4.96E-003 




2.36 


7.09E-003 


1.98 


5.57E-003 


2.63 




64 


1.16E-003 




2.09 


1.90E-003 


1.90 


1.18E-003 


2.24 




128 


2.76E-004 




2.08 


4.25E-004 


2.16 


3.24E-004 


1.86 




256 


6.73E-005 




2.04 


9.32E-005 


2.19 


9.67E-005 


1.75 


8W 


16 


2.60E-002 






2.81E-002 




3.37E-002 






32 


5.19E-003 




2.32 


7.28E-003 


1.95 


5.44E-003 


2.63 




64 


1.22E-003 




2.09 


1.88E-003 


1.95 


1.37E-003 


1.99 




128 


2.96E-004 




2.05 


4.02E-004 


2.22 


3.45E-004 


1.99 




256 


7.29E-005 




2.02 


8.40E-005 


2.26 


8.79E-005 


1.97 



We begin by considering two dimensional propagation choosing tana = 2, (3 = in accordance with [18L 
Il5l.[l7j]. Computations are carried out for exactly one wave period (t = T = l/y/E) on the computational box 
[0, 1] x [0, 1/2] with N x x N y points, where N y = N x /2. Errors, computed as ^ei(B x ) 2 + e\{B y ) 2 + ei{B z ) 2 , 
are reported in Tablc[T]and plotted as function of the mesh size, N x = 16, 256, in the left panel of Fig[TJ 

Selected schemes (CT, GLM and 8W) produce comparable errors and show essentially second-order 
accuracy. We notice that decreasing the Courant number to C a = 0.4 has the effect of slightly reducing the 
errors for GLM at large resolution but not for CT and 8W. From Table [TJ in fact, one can see that, when 
N x = 256, the error is reduced from ~ 8.5 • 10~ 5 to ~ 7 • 10 -5 for GLM, while it grows from 6.7 • 10 -5 to 
9.3 • 10~ 5 for the CT scheme. 
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We have found that the solution is very weakly dependent on the a parameter and the errors are 
minimized when a = 0. Besides, we repeated the computations with the EGLM formulation and observed 
essentially the same level of accuracy with no particular improvement over GLM. 

4-1-2. Three- Dimensional Propagation 

In three dimensions we follow (ljj and set tan a = tan f3 = 2 so that the resulting computational box is 
given by a; £ [0, 1], y, z € [0, 1/2] discretized on N x x N x /2 x N x /2 grid points. Computations are followed 
for one wave period (T = 1/3) and repeated, with C a = 0.4, on increasingly finer grids corresponding to 
N x = 16,32,64,128,256. The right panel in Figure [1] shows that all schemes meet the expected order of 
accuracy providing comparable errors, as found in the two-dimensional case. A more quantitative comparison 
can be made by inspecting the last two columns of Table [TJ where one can see that GLM performs slightly 
better than the other schemes. 



4-2. Nonlinear smooth flow 

In the next example we consider the evolution of a fully nonlinear smooth flow where, unlike the previous 
example, all waves (linear and nonlinear) are triggered. Following p?! [l[ we specify a periodic computational 
box in Cartesian coordinates, spanning from —1 to 1 in the x and y directions with initial conditions given 
by 



1 



1 



(V X , Vy) 



(B X , By) 



■ + - sin(Tre) + - cos(?ry) , 
1 + — sin(7ry) + — cos(7rx), 1 + — sin(7ra;) + — cos(7r?/) 



(38) 



where p = 1/4, while v z = B z = 0. Integration terminates at t = 0.2, before the formation of any 
discontinuous feature. A resolution study is carried out for all schemes and compared to a reference solution 
obtained on 2048 2 zones with the CT scheme. The error, shown in Fig. ^ as a function of the number of 
cells, is computed as a quadratic mean of the L\ norm errors (given by Eq. I32p of the primitive variables. All 
schemes are second-order accurate with comparable errors, with the GLM approach giving slightly better 
results than the others at the largest resolution (256 zones). 



4-3. Shock Tube Problems 

One dimensional shock tubes have proven to be valuable benchmarks in order to assess the ability of the 
scheme to capture both continuous and discontinuous flow features. The rotated multidimensional versions 
considered in the following may be used to check the strength of the numerical method in preserving the 
original planar symmetry through an oblique propagation. 

4-3.1. Two-dimensional shock tube 

In the first shock tube, taken from [26j . we consider an initial discontinuity with left and right states given 
by (p, v 1 ,v 2 ,B 1 ,B 2 ,p) L = (1, 10, 0, 5/V4^, 5/^4^, 20) and (p, Vl ,v 2 , B u B 2 ,p) R = (1, -10, 0, 5/V^, 5/^4^, 1) 
respectively. The subscripts "1" and "2" give the directions perpendicular and parallel to the initial discon- 
tinuity. The initial condition is then rotated on a Cartesian grid (x, y) using the transformation defined by 
Eq. with a = tan" 1 2 and ft = 7 = 0. 

Since the magnetic field is initially uniform, V • B = is trivially ensured at t = 0. The computational 
domain spans from to 1 in the x direction and from to 2/N x in the y direction with N x x 2 computational 
zones. Outflow boundaries are set at the rightmost and leftmost sides of the box whereas for any flow 
variables q at the upper and lower boundaries we impose the translational invariance q{i,j) = q(i±5i,j±6j) 
where (Si,Sj) = (2, —1) with the plus (minus) sign holding at the upper (lower) boundary. Computations 
terminate before the fast shocks reach the boundaries, at t = 0.08 cos a. 
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Table 2: One dimensional L\ (xlO 2 ) norm error for the two-dimensional shock tube. 



P Vi V 2 B x B 2 p 

8W 2.7 8.6 1.9 9.6 6.2 94.5 

CT 2.6 8.5 1.5 0.4 4.7 93.0 

GLM 2.6 8.4 1.4 0.4 4.3 90.5 

EGLM 3.2 8.3 1.3 0.4 5.1 96.4 



Fig [3] shows the primitive variable profiles obtained with the conservative GLM-MHD scheme. The 
resulting wave pattern is comprised of two outermost fast shocks (at x\ ~ 0.12 and x\ ~ 086) enclosing two 
slow magneto-sonic waves and a contact mode at x\ ~ 0.56. We see that all discontinuities are captured 
correctly although some spurious oscillations are visible in the transverse velocity profile in proximity of the 
fast shocks. Similar features are also evident in the paper by Toth (2^| and with the CT scheme (not shown 
here) . 

We have repeated the same test with the four different schemes described at the beginning of this section 
and compared the results against a one-dimensional reference solution obtained at higher resolution (1024 
cells) up to t — 0.08. Tabic [2] gives the errors, using the one-dimensional L\ norm, of the primitive variables 
for the 8W, CT, GLM and EGLM schemes. While errors in density, velocity and pressure arc very similar 
for all schemes, the longitudinal component of the magnetic field [B\] shows substantially large deviations 
with the 8W scheme. This is further illustrated in Fig 2] where, in accordance with (26[, we find that the 
8-wave formulation results in erroneous jump conditions in the normal component of the field. On the other 
hand, both the GLM and the non conservative EGLM schemes behave as well as CT on this particular test 
without producing spurious jump conditions. 

Finally, in the left panel in Fig. [7] we plot, as a function of a, the L\ norm errors in B\ at different 
resolutions, N x = 128, 256, 512, for both the GLM (black) and EGLM (red) formulations. The plots show 
a weak dependence on the a parameter and errors are minimized for a w 0.5, independently of the mesh 
resolution, for both schemes. Also, owing to the presence of shock waves, the order of convergence is 
approximately one. 

4-3.2. Three-dimensional shock tube 

For the three dimensional version we follow [lj| and set the initial left and right states to 

2 3 6 2 \ T 
V L = I 1.08, 1.2,0.01,0.5,==, ====,0.95 for xi<0, 

V4^ V4^V4^ J 

2 4 2 



V fl = 1,0,0,0,-=,-=,-=,! for Xl >0, 

a/47T V47T a/47T 



wherc V = (p, v\, v 2 , i>3, B\, B 2 , B%,p) is the vector of primitive variables. The coordinate transformation 
used for the 3D rotation is given by Eq. ([33)) where the rotation angles a and (3 are chosen in such a 
way that an integer shift of cells satisfies, for any flow quantities q, the translational invariance expressed 
by g(x + s) = g(x), where s is a Cartesian vector orthogonal to xi and thus xi(x + s) = Xi(x). This 
condition follows from the fact that the solution is a function of x\ alone and thus invariant for translations 
transverse to this direction, providing a convenient way to assign boundary conditions in the (x, y, z) system 
of coordinates. By choosing tana = —ri/r 2 and tan/3 = ri/r$ together with s = (n x Ax, n y Ay, n z Az), one 
can show that the three shift integers n x ,n y , n z must obey 

7*1 7*1 

n x — n y \-n z — =0, (40) 

r 2 r 3 

where Ace = Ay = Az has been assumed and (r*i, r2, 7*3) = (1,2,4) will be used. The computational domain 
consists of [768 x 8 x 8] zones and spans (—0.75, 0.75) in the x direction while y, z £ [0, 0.015625]. 
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Tabic 3: Li (xlCT 4 ) error for the 3D Shock Tube 



P Vl V2 V 3 B 1 B 2 B 3 p 

8W 3.0 2.0 4.8 4.7 3.6 4.5 5.1 5.0 

CT 3.1 2.4 4.2 4.4 0.5 5.3 5.4 5.5 

GLM 2.9 2.3 3.6 4.3 0.5 4.7 5.4 5.1 

EGLM 3.5 2.5 4.3 4.8 0.5 5.3 5.9 7.3 



In Fig. [5] we plot the primitive variable profiles for the GLM scheme at t = 0.02 cos a cos 7. In accordance 
with the one dimensional solution (see also [Io|). we observe the formation of a structure involving a contact 
discontinuity separating two fast shocks, two slow shocks and a pair of rotational discontinuities. The 
three-dimensional integration reproduces the correct behavior of all waves and the error in the longitudinal 
component of the field (Bi in Fig exhibits small spurious oscillations about the same order of the CT 
scheme (see also, for instance, Fig. 7 in [lil]). 

A quantitative estimate of the error (using the one-dimensional L\ norm error) is obtained by comparing 
the three-dimensional results with a one-dimensional reference solution computed on 1024 zones until t = 
0.02. The comparison, extended to the four selected integration schemes, is given in Table [3] We notice 
that the CT, GLM and EGLM schemes all yield errors of the same order of magnitude (typically 10 -4 ). 
Beware that these computations may be susceptible to small variations depending on implementation details 
(e.g. limitcr, Courant number, etc.) and thus give a representative estimate of the error. For instance, the 
implementation of the CTU-CT scheme in the PLUTO code [2(| is similar, although not exactly equivalent, 
to that of fl6j who instead use piecewise parabolic reconstruction. Nevertheless, we have ascertained that 
the 8W scheme always performs the worst and the discrepancy becomes particular evident by looking at 
the longitudinal component of the field where the 8W scheme yields, once again, incorrect (although smaller 
than the previous 2D case) jumps. This is better illustrated in Fig. [6l where we compare the profiles of B\ 
for the four selected numerical schemes. We stress that, despite its non-conservative character, the EGLM 
formulation does not seem to produce incorrect jump conditions or wrong shock propagation speeds. 

A resolution study, shown in the right panel of Fig [71 demonstrates that errors produced by the GLM 
and EGLM formulations are very much comparable and only weakly dependent on the a parameter. Both 
schemes report a minimum at a m 0.005 — 0.01 regardless of the resolution, and the inferred order of 
convergence is approximately one as expected for solutions involving shock waves. 

4-4- Magnetic Field Loop Advection 

This problem consists of a weak magnetic field loop being advected in a uniform velocity field. Since 
the total pressure is dominated by the thermal contribution, the magnetic field is essentially transported as 
passive scalar. 



4-4-1- Two-dimensional advection 

Following 15|, 14,[l7|, we employ a periodic computational box defined by x G [—1, 1] and y 6 [—0.5, 0.5] 



discretized on N x x N x /2 grid cells (N x = 128). Density and pressure are initially constant and equal to 1. 
The velocity of the flow is given by v — (Vo cos a, Vq sin a, 1) with Vb = y/5, sin a = l/y/E and cos a = 2/Vo. 
The magnetic field is defined through its magnetic vector potential as 

( Ao(R-r) if r<R, 
A z = \ ~ (41) 

if r > R , 



where Ao = 10~ 3 , R = 0.3 and r = \Jx 2 + y 2 . The simulations are allowed to evolve until t = 2 ensuring 
the crossing of the loop twice through the periodic boundaries. 

In Fig. [5] we show the magnetic energy density for the 8W, GLM and CT schemes using C a = 0.8 (top) 
and C a = 0.4 (bottom), along with the field lines shape. The circular shape of the loop is best preserved 
with the CT and GLM schemes while some distortions are visible using the 8 wave formulation. Using 
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C a = 0.4 with the GLM scheme yields slightly better results, while the CT docs not seem to be affected by 
the choice of the Courant number. 

The time- history of the magnetic energy density (left panel in Fig[5]) reveals that the numerical dissipation 
is essentially similar for all schemes, being smaller at larger Courant numbers. At the quantitative level, our 
results are similar and in good agreement with those of other investigators (e.g., .[l5, 14. Il7|). 

The ability of the GLM scheme in preserving the divergence-free condition is monitored by checking the 
growth of B z in time: owing to a non-vanishing z component of velocity, in fact, we expect B z to grow in 
time with a rate oc w z V • B as seen from the induction equation. In the middle panel of Fig. [5] we plot the 
volume-averaged value of B z as a function of time for N x = 64, 128, 256. The nominal value is ~ 10~ 3 of the 
initial field strength, decreasing with resolution. Notice that the observed order of convergence is ~ 0.6 — 0.7 
and thus sub-linear as expected for a linearly degenerate wave in Godunov-type schemes, in accordance with 
the results of [f|. 

Computations carried with different values of a reveal that divergence errors are minimized for a > 0.01 
while errors in B z become smallest for a ss 0.01 (right panel in Fig [5]). Despite this may generate some 
ambiguities in prescribing an optimal a value, however, we see that its choice does not significantly affect 
the error and thus constitutes a minor effect on the solution. 

4-5. Three-dimensional field loop advection 

A three-dimensional extension can be obtained by rotating the previous 2D configuration around one axis 
using the coordinate transformation given by Eq. ([55]) with a = and 7 = tan -1 1/2, see [lq ]. Even though 
the loop is rotated only around one axis, the velocity profile (v x , v y , v z ) = (1,1,2) makes the test intrinsically 
three-dimensional. We consider the computational box —0.5 < x < 0.5, —0.5 < y < 0.5, —1.0 < z < 1.0, 
resolved on a N x N x 27V grid. Boundary conditions are periodic in all directions. 

A three-dimensional rendering of the magnetic energy density is shown in Fig. [TU]for the selected schemes 
while relevant quantities are plotted in the three panels of Fig 111! All schemes show a similar amount of 
numerical dissipation, in agreement with the results of [l6j ] . 

As for the 2D case, it is useful to check the growth of the magnetic field component B% = (— B X + 2B Z ) / \/Z 
orthogonal to the original (3^,2:2) plane where the loop is two-dimensional. Analytically, the magnetic field 
component in this direction is a trivial constant of motion since 

m + -0. (-) 



dt \ dx\ 0x2 

The numerical integration in the rotated (x, y 1 z) Cartesian frame, however, preserves this condition only to 
some accuracy which strongly reflects the ability of the scheme in controlling the divergence-free constraint 
(this is true for all presented numerical methods). The middle panel in Fig I 111 shows the volume- integrated 
value of I-B3I, normalized to the initial field strength Bq = 10~ 3 for three different resolutions TV = 32, 64, 128. 
Our results reveal that the value of B3 grows slowly in time while remaining reasonably small. The conver- 
gence rate (« 0.6 — 0.7) is approximately the same as the one observed in the 2D case. 

The dependency on a is illustrated in the right panel Fig [11] showing that divergence errors are progres- 
sively reduced for a > 0.03 although this has very little effect on the growth of B3. 

4-6. Two-dimensional Rotor problem 

The rotor problem consists of a dense disk rotating in a static medium threaded by an initially uniform 
magnetic field. As the rotor spins, the magnetic field gets wrapped around the disk creating torsional Alfvcn 
waves, stemming from the rotating disk and moving towards the surrounding gas. This interaction slows 
down the disk by extracting angular momentum. On the other hand, the build-up of magnetic pressure 
around the rotor causes its compression. 
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Table 4: Parameter sets used for the first and second versions of the three-dimensional blast wave problem. 





Pin 


Pout 


B 


9 


^0 


^stop 


Test 1 


10 2 


1 


10 


tt/4 


0.125 


0.02 


Test 2 


10 4 


1 


100 





0.1 


2.5- 10~ 3 



We initialized the problem on the Cartesian box x,y £ [— \, \\ with outflow boundary conditions and 
use 400 2 grid points. The primitive variable profiles at the beginning of the simulation are given by 

(10, -cjy, ujx) if r<r , 

(p,v m ,v y ) = i (l + 9f,-f(jy^-,fujx^j ifr <r<n, (43) 

(1,0,0) if r>r x , 



where ui = 20, ro = 0.1, r\ = 0.115, r = y/ x 2 + y 2 and the taper function is / = (ri — r)/(ri — ro). 
Thermal pressure is initially uniform and equal to one (T = 1.4 is used). The magnetic field has only one 
non- vanishing component, B x = 5/\/47r. 

The maps of density, magnetic energy and sonic Mach number are displayed in Fig. [l_2]at t = 0.15 for the 
GLM and the CT schemes, when the torsional Alfen waves have almost reached the outer boundaries. The 
strength of the scheme is also measured by its ability to preserve the circular shape of the sonic Mach number 



profile in the central region, an essential feature of the solution, 17[. This is better shown in Fig. [T31 where 
an enlargement of the central region reveals that the GLM and CT schemes have developed extremely similar 
Mach number contours and the absence of spurious peaks (that would be caused by pressure undershoots) 
advocates towards the validity of the scheme. 

4-7. Three Dimensional Blast Wave 

The MHD blast wave problem has been specifically designed to show the scheme ability to handle 
strong shock waves propagating in highly magnetized environments, see for instance 1, H [H, U3. 



Depending on the strength of the magnetic field, it can become a rather arduous test leading to unphysical 
densities or pressures if the divergence-free condition is not properly controlled and the scheme does not 
introduce adequate dissipation across oblique discontinuous features. Here, we consider a three-dimensional 
configuration on the unit cube [—1/2, 1/2] 3 discretized on 200 3 computational zones. The medium is initially 
at rest (v = 0) and threaded by a constant uniform magnetic field lying in the xz plane and forming an 
angle 9 with the vertical z direction, B = Bo (sm9x + cos6z). A spherical region of high thermal pressure 
is initialized, 

{Pin for J x 2 + y 2 + z 2 < r , 

(44) 
p ou t otherwise . 

We consider two different versions of the same test problem with parameters given in Table |H In the first 
one, taken from [l6j |. the field forms an angle 9 = 7r/4 with the z axis and the largest magnetization achieved 
outside the sphere is (3 = 2p out /B 2 = 2 • 10 -2 . In the second version, we follow [29[ and adopt a a larger 
field strength (with = 0) yielding a more severe configuration with = 2- 10 -4 . 

The over-pressurized spherical region sets a blast wave delimited by an outer fast forward shock prop- 
agating (nearly) radially, sec Fig [TJ] and [TBI Magnetic field lines pile up behind the shock in the direction 
transverse to the initial field orientation (8 = 7r/4 and 9 — for the two cases) thus building a region of higher 
magnetic pressure. In these regions the shock becomes magnetically dominated and only weakly compressive 
(Sp/p ~ 1.2 in both cases). The inner structure is delimited by an oval-shaped slow shock adjacent to a 
contact discontinuity and the two fronts tend to blend together as the propagation becomes perpendicular 
to the field lines. The magnetic energy increases behind the fast shock and decreases downstream of the 
slow shock. The resulting explosion becomes highly anisotropic and magnetically confined. 
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Computed results for the first configuration are shown in Fig 1141 where we display linearly scaled maps 
of gas pressure, magnetic and kinetic energy densities for the GLM scheme (top), EGLM (middle) and 
CT schemes (bottom). The computations are in excellent agreement and no noticeable difference can be 
discerned from the images. Moreover, our results favorably compare to those of [T^J. To further ascer- 
tain the validity of the non-conservative EGLM scheme, we plot, in Fig [T5j one-dimensional slices (along 
the x direction in the yz mid-plane) showing the density and pressure obtained with the EGLM and CT 
integrations. 

Computations for the second configuration could be obtained only with the EGLM scheme, since the CT 
scheme failed even with a minmod limiter (/3 = 1 in Eq. I18p . In Fig. [16] we plot contour levels for density, 
pressure, velocity and magnetic energy. These results comply with those of [29[ who used a CT scheme 
together with a Runge-Kutta time stepping and an HLL Ricmann solver. They also share similarities with 
the 2D strong field case discussed in [17| who used a different implementation of the CT scheme. Partially 
owing also to the increased resolution (200 3 instead of 144 3 ) our CTU-EGLM algorithm shows considerably 
reduced numerical diffusion while being robust in keeping sharp profiles of the discontinuities. 

5. Conclusions 

A second-order, cell-centered numerical scheme for the solution of the MHD equations in two and three 
dimensions has been proposed. Fully unsplit integration resorts to the Corner Transport Method of Colella 
and the divergence-free condition is controlled by using a constrained formulation of the MHD equations 
where the induction equation is coupled to a generalized Lagrange multiplier (GLM, [12]). The system 
is hyperbolic, easy to implement and does not require expensive cleaning projection steps associated with 
the solution of elliptic problems. The GLM scheme is fully conservative in mass, momentum, energy and 
magnetic induction, although we have also considered a slightly modified variant (EGLM) which infringes 
momentum and energy conservation. 

In order to assess the reliability and accuracy of the schemes we have performed a number of code 
benchmarks on standard two- and three-dimensional MHD test problems. Results have been compared with 
two different numerical schemes: a non-conservative cell-centered method based on the 8-wave formulation 
(8W, [22j|) and the constrained transport (CT) method where the magnetic field has a staggered collocation. 
Both the GLM and EGLM schemes give excellent results in terms of accuracy and robustness and do not 
show, in the tests presented here, any evidence for incorrect jump conditions or wrong wave propagation, 
as found for the eight wave formulation (in agreement with Toth [26j). This has been verified on problems 
involving discontinuous waves and holds true for both the conservative GLM formulation and the EGLM 
variant which breaks momentum and energy conservation. In this perspective, our results seem to indicate 
that the presence of source terms in the equations does not necessarily lead to erroneous jumps. Instead, 
we have found the non-conservative formulation to be more robust for problems involving the propagation 
of oblique strongly magnetized shocks. Although, this behavior may be attributed to discretization, such a 
study is beyond the scope of the present paper. The comparison has also revealed an excellent quantitative 
agreement with the CTU-CT scheme (in the version 

of HEE3) 

showing errors with comparable magnitude 
and similar order of convergence while retaining the desired robustness and stability. 

For these reasons, we believe that the proposed CTU-GLM and CTU-EGLM schemes provide excellent 
competitive alternatives to modern staggered-mesh algorithms while being considerably easier and more flex- 
ible in their implementations. Owing to the cell-centered collocation of all of the flow fields, the CTU-GLM 
scheme can be easily generalized to resistive MHD, adaptive and/or unstructured grids and to higher than 
second-order spatially-accurate numerical schemes. Some of these issues will be presented in forthcoming 
papers. 

A. Characteristic Decomposition of the GLM-MHD Equations 

The 9x9 matrix A x of the primitive MHD equations introduced in ^3.1.1l can be decomposed as A x = RAL 
where A = diag(A fc ) contains the eigenvalues (see Eq. [T4"]) while the rows of L and columns of R are the 
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corresponding left and right eigenvectors of A x , respectively. Adopting the scaling of [22| we define 



a 



dj-c 2 



c 2 ~c 2 



and 



0v 



B„ 



B, 



B 2 +B 2 



B 2 + B 2 

y z 



(45) 
(46) 



where a = y/Tp/p denotes the speed of sound. With this notation, the right eigenvectors in matrix form 
will be given by 
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where S = sign(B x ). On the other hand, the left eigenvectors are 
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10 100 10 100 

N N 

Figure 1: L\ norm errors for the 2D (left) and 3D (right) circularly polarized Alfven wave test problem. Each symbol refers 
to results obtained with the GLM (plus sign), CT (square) and Powell's eight wave (rhombus) methods, while the dotted line 
gives the ideal second-order convergence slope. The Courant number C a = 0.8 and the final time step is l/v5 (left) and 1/3 
right. 
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Figure 2: L\ norm errors for the non-linear, smooth flow test problem at t = 0.2. The different symbols refer to computations 
carried out with the GLM (plus signs), EGLM (ex signs), CT (squares) and Powell's eight wave method (rhombus) with 
Courant number C a = 0.8. The dotted line gives the ideal second-order convergence slope. 
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Figure 3: Primitive variable profiles for the 2D shock tube problem at t = 0.08 cos a, along the rotated direction xi- The 
symbols correspond to the CTU-GLM solution whereas the solid lines represent the reference solution. From top to bottom 
and left to right, density, thermal pressure, velocity components and magnetic field components (parallel and perpendicular 
with respect to the xi direction) are displayed. 
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Figure 4: The parallel magnetic field component for the four schemes. Concordantly with the results of 26] the 8 wave formalism 
fails to capture the correct jumps. This problem is absent in the results of the other schemes and the field component remains 
close to the expected value 5/%/~4tt away from discontinuities. Spikes are found in proximity of shock waves and are of the same 
order of magnitude for GLM, EGLM and CT schemes. 
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Figure 5: Primitive variable profiles for the 3D shock tube problem at t = 0.02 cos a cos 7, along the rotated direction xi 
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Figure 6: Comparison of the parallel component of the magnetic field for the 3D shock tube test. As in the 2D case, the error 
is minimal for all schemes with the exception of the 8-wave formalism. The latter fails to capture correctly the jump but the 
error is less prominent than the 2D case. 
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Figure 7: L\ norm errors of B\ (the magnetic field component in the direction orthogonal to the initial discontinuity) as functions 
of a = Ahch/Cp for the 2D (left) and 3D (right) shock tube problem. The different symbols correspond to computations 
carried at different mesh resolutions: N x = 128, 256, 512 (in 2D) and N x = 384, 768, 1536 in 3D. Black and red symbols refer 
to results obtained with the GLM and EGLM formulations, respectively. The Courant numbers were 0.8 and 0.4 for 2 and 3D 
computations, respectively. 
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Figure 8: From left to right: magnetic energy density for the 2D field loop problem at t = 2 for the 8W, GLM and CT schemes. 
Results have been computed with CFL numbers of 0.8 (top) and 0.4 (bottom). Overplotted are 9 isocontours of A z , between 
10" 5 and 10" 3 . 



<B 2 /B 2 >(t) <IB Z I/B >(t) Optimal a 
1.001 J I I I 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.00 0.02 0.04 0.06 0.08 0.10 
t t a 



Figure 9: Leftmost panel: time evolution of the volume-integrated magnetic energy density (normalized to its initial value) for 
the 2D field loop advection problem. The black and red lines correspond, respectively, to computations carried with C a = 0.4 
and C a = 0.8. Middle panel: volume-averaged value of \B Z | (normalized to the initial value Bo = 10 — 3 ) as a function of time for 
three different grid resolutions (256, 128 and 64 corresponding to stars, "x" and plus signs). Rightmost panel: volume-averaged 
values of | V • B | and | B z | for different values of the a parameter controlling monopole damping at the resolution N x = 128 
points. 
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Figure 10: Magnetic energy density for the 3D field loop problem at t = 1 at the resolution of 128 X 128 X 256. From left to 
right: results obtained with the 8W, GLM and CT schemes. 
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Figure 11: Same as Fig [9] for the 3D field loop advection test. From left to right: time history of the (normalized) volume- 
integrated magnetic field energy, (normalized) average value of \B^\ (magnetic field component orthogonal to the original 2D 
plane) and volume averages of |V ■ B3I and \B%\ as functions of the a parameter. 
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Figure 13: A zoom in the central region of the rotor problem at time t = 0.15, showing 20 levels (0 < M < 4) of contour 
profiles of the sonic Mach number. Results for the GLM and CT schemes are shown on the left and right panels, respectively. 
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Figure 14: Two dimensional cuts in the xz plane of gas pressure, magnetic and kinetic energy densities for the GLM (top), 
EGLM (middle) and CT (bottom) schemes, at t = 0.02 for the first blast wave problem. Pressure values range from 1.0 (white) 
to 42.4 (black). The magnetic energy ranges from 25.2 (white) to 64.9 (black) while the kinetic energy density spans from 0.0 
(white) to 33.1 (black). 



29 




-0.4 



-0.2 0.0 



0.2 



10 - 



0.4 




Figure 15: Density and pressure profiles, the latter on logarithmic scale, along x at y, z = 0, at time t = 0.02. Results obtained 
with the CT and EGLM schemes are shown using box and cross symbols, respectively. 
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Figure 16: Density, pressure, velocity and magnetic energy contours (30 levels) for the CTU-EGLM scheme at t = 2.5 ■ 10 — 3 
in the xz plane. Density values range from 0.18 to 3.2 while pressure spans from 0.9 to 2290. The absolute value of velocity 
ranges from 0.0 to 47 while the magnetic energy spans from 2817 to 5932. 
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